#ifndef EVECTOROPS2_HPP_
#define EVECTOROPS2_HPP_

/* The following code example is taken from the book
 * "C++ Templates - The Complete Guide"
 * by David Vandevoorde and Nicolai M. Josuttis, Addison-Wesley, 2002
 *
 * (C) Copyright David Vandevoorde and Nicolai M. Josuttis 2002.
 * Permission to copy, use, modify, sell and distribute this software
 * is granted provided this copyright notice appears in all copies.
 * This software is provided "as is" without express or implied
 * warranty, and with no claim as to its suitability for any purpose.
 */
/* modified by ofgr 21.03.2008 */
#include <math.h>
#include "usertypes.hpp"

//*****************
// unary operators
//*****************

// keeping sign
template <typename T, typename R2>
EVector<T,A_Plus<T,R2> >
operator+ (EVector<T,R2> const& b) 
{
  return EVector<T, A_Plus<T,R2> > (A_Plus<T,R2>(b.rep()));
}

// changing sign
template <typename T, typename R2>
EVector<T,A_Minus<T,R2> >
operator- (EVector<T,R2> const& b) 
{
  return EVector<T, A_Minus<T,R2> > (A_Minus<T,R2>(b.rep()));
}



//**********
// ADDITION
//**********
// addition of two EVectors
template <typename T, typename R1, typename R2>
EVector<T,A_Add<T,R1,R2> >
operator+ (EVector<T,R1> const& a, EVector<T,R2> const& b) 
{
  return EVector<T,A_Add<T,R1,R2> >
    (A_Add<T,R1,R2>(a.rep(),b.rep()));
}

// addition of scalar and EVector
template <typename T, typename R2>
EVector<T,A_Add<T,A_Scalar<T>,R2> >
operator+ (T const& s, EVector<T,R2> const& b) 
{
  return EVector<T, A_Add<T,A_Scalar<T>,R2> >
    (A_Add<T,A_Scalar<T>,R2>(A_Scalar<T>(s),b.rep()));
}

// addition of EVector and scalar
template <typename T, typename R1>
EVector<T, A_Add<T,R1,A_Scalar<T> > >
operator+ (EVector<T,R1> const& a, T const& s) 
{
  return EVector<T,A_Add<T,R1,A_Scalar<T> > >
    (A_Add<T,R1,A_Scalar<T> > (a.rep(), A_Scalar<T>(s)));
}



//*************
// SUBTRACTION
//*************
// subtraction of two EVectors
template <typename T, typename R1, typename R2>
EVector<T,A_Sub<T,R1,R2> >
operator- (EVector<T,R1> const& a, EVector<T,R2> const& b) 
{
  return EVector<T,A_Sub<T,R1,R2> >
    (A_Sub<T,R1,R2>(a.rep(),b.rep()));
}

// subtraction of scalar and EVector
template <typename T, typename R2>
EVector<T,A_Sub<T,A_Scalar<T>,R2> >
operator- (T const& s, EVector<T,R2> const& b) 
{
  return EVector<T, A_Sub<T,A_Scalar<T>,R2> >
    (A_Sub<T,A_Scalar<T>,R2>(A_Scalar<T>(s),b.rep()));
}

// subtraction of EVector and scalar
template <typename T, typename R1>
EVector<T, A_Sub<T,R1,A_Scalar<T> > >
operator- (EVector<T,R1> const& a, T const& s) 
{
  return EVector<T,A_Sub<T,R1,A_Scalar<T> > >
    (A_Sub<T,R1,A_Scalar<T> > (a.rep(), A_Scalar<T>(s)));
}

//****************
// MULTIPLICATION
//****************
// multiplication of two EVectors
template <typename T, typename R1, typename R2>
EVector<T, A_Mult<T,R1,R2> >
operator* (EVector<T,R1> const& a, EVector<T,R2> const& b) 
{
  return EVector<T,A_Mult<T,R1,R2> >
    (A_Mult<T,R1,R2>(a.rep(), b.rep()));
}

// multiplication of scalar and EVector
template <typename T, typename R2>
EVector<T, A_Mult<T,A_Scalar<T>,R2> >
operator* (T const& s, EVector<T,R2> const& b) 
{
  return EVector<T,A_Mult<T,A_Scalar<T>,R2> >
    (A_Mult<T,A_Scalar<T>,R2>(A_Scalar<T>(s), b.rep()));
}

// multiplication of EVector and scalar
template <typename T, typename R1>
EVector<T, A_Mult<T,R1,A_Scalar<T> > >
operator* (EVector<T,R1> const& a, T const& s) 
{
  return EVector<T,A_Mult<T,R1,A_Scalar<T> > >
    (A_Mult<T,R1,A_Scalar<T> > (a.rep(), A_Scalar<T>(s)));
}


//****************
// SCALARPRODUCT
//****************
// scalarproduct of two EVectors
template <typename T, typename R1, typename R2>
T inner_product (EVector<T,R1> const& a, EVector<T,R2> const& b) 
{
  //assert (a.size () == b.size ());
  T res = T();

  for (myindex i=0; i<a.size (); ++i) 
    res += a[i]*b[i];
  return res;
}

//****************
// ABS
//****************
// absolute value / norm of an EVectors 
template <typename T, typename R1>
T abs (EVector<T,R1> const& a)
{
  return sqrt (inner_product(a,a));
}

//****************
// CROSSPRODUCT
//****************
// evaluation must not be realized by expression template 
// due to row-wise execution (which can give wrong results)
template <typename T, typename R1, typename R2>
EVector<T> 
crossproduct (EVector<T,R1> const& x, EVector<T,R2> const& y) 
{ 
  assert (x.size () == 3 && y.size () == 3);
  EVector<T> z(3);
  z[0] = x[1] * y[2] - x[2] * y[1];
  z[1] = x[2] * y[0] - x[0] * y[2];
  z[2] = x[0] * y[1] - x[1] * y[0];

  return z;
}

#endif //EVECTOROPS2_HPP_
